Time evolution of a quantum many-body system: transition from integrability to 

ergodicity in thermodynamic limit 
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Numerical evidence is given for non-ergodic (non-mixing) 
behavior, exhibiting ideal transport, of a simple non- 
integrable many-body quantum system in the thermodynamic 
limit, namely kicked t—V model of spinless fermions on a ring. 
However, for sufficiently large kick parameters t and V we re- 
cover quantum ergodicity, and normal transport, which can 
be described by random matrix theory. 

PACS numbers: 05.30.Fk, 05.45.+b, 72.10.Bg 



A simple question is addressed here: "Do mixed quan- 
tum many-body systems, which are neither integrable 
nor ergodic, exist in the thermodynamic limit?" While 
it is clear that integrable systems are rather exceptional 
it is an important open question whether a finite generic 
perturbation of an integrable system becomes ergodic or 
not in the thermodynamic limit (TL), size — » oo and 
fixed density. It is known that local statistical prop- 
erties of quantum systems with few degrees of freedom 
whose classical limit is completely chaotic — ergodic are 
universally described by random matrix theory, while in 
the other extreme case of integrable systems Poissonian 
statistics may typically be applied jj],^]. This statement 
has recently been verified numerically also for integrable 
and strongly non-integrable many-body systems of inter- 
acting fermions || which do not have classical limit. 

Having lost the reference to classical dynamics we re- 
sort to the definition of quantum ergodicity (also termed 
quantum mixing) 0] as the decay of time-correlations 
(A(t)B(O)) — (A)(B) of any pair of quantum observables 
A and B in TL, taking the time-limit r — > oo in the 
end. In B many-body system of interacting bosons has 
been studied and it has been shown that quantum ergod- 
icity corresponds to strongly chaotic (classically ergodic) 
dynamics of associated non-linear mean-field equations. 
As a consequence of linear response theory, quantum er- 
godicity also implies normal transport and finite trans- 
port coefficients (such as dc electrical conductivity). On 
the other hand, integrable systems, which are solvable by 
Bethe ansatz or quantum inverse scattering, are charac- 
terized by (infinitely many) conservation laws and are 
thus non-ergodic. It has been pointed out recently |j| 
that integrability implies nonvanishing stiffness, i.e. ideal 
conductance with infinite transport coefficients (or ideal 
insulating state). As we argue below, any deviation 
from quantum ergodicity generically implies nonvanish- 
ing long-time current auto-correlation and therefore infi- 
nite transport coefficient. Since generic non-integrable 



systems of finite size (number of degrees of freedom) 
are non-ergodic (obeying mixed statistics smoothly in- 
terpolating from Poissonian to random matrix results) 
it is thus an important question if and when such non- 
ergodicity can survive TL. 

In this Letter we introduce a family of simple many- 
body systems smoothly interpolating between integrable 
and ergodic regime, namely kicked t- V model (KtV) 
of spinless fermions with periodically switched nearest 
neighbor-interaction on a 1-dim lattice of size L and pe- 
riodic b.c. L = 0, with time-dependent hamiltonian 

L-l 

H ^ = E [-H c h+ 1 + h - c -) + S p( T ) Vn J n J+i] > (!) 

and give numerical evidence for the existence of mixed 
regime (reader should not confuse it with mixing) in TL 
by direct simulation of the time evolution. Cj,Cj,n,j are 
fermionic creation, annihilation and number operators, 
respectively, and S p (r) — X)m=-oo < H r ~ m )- Deviations 
from quantum ergodicity (or mixing) are characterized 
by several different quantities as described below. 

KtV model ([!]) is a many-body analogue of popular 1- 
dim non-integrable kicked systems Q such as e.g. kicked 
rotor: its evoulution (Floquet) operator over one period 

U = Texp(— i J Q+ dTH(r)) (h = 1), factorizes into the 
product of kinetic and potential part 



U = cxp 



-iV 



L-l 

E 

3=0 



rijnj+i exp 




cos(sfc 



(2) 



where s = 2-k/L. Flux parameter 
to introduce a current operator J = 



cf> is used in order 
= {%/t)Wd^U\^a = 
X)fc=o shx(sk)hki elsewhere we put <f> := 0. Tilde de- 
notes the operators which refer to momentum variable k, 
Ck = L~ 1 / 2 J2f=o exp(isjk)cj, h k = c\c k . KtV model is 
integrable if either t = 0, or V = (mod 2tt), or tV — > 
and t/V finite (continuous time t-V model, see e.g. ||), 
while for t ~ V ~ 1 it is expected to be non-integrable, 
either quantum ergodic or mixed. 

We expect that unitary many-body quantum maps, 
such as (||), also mimick the dynamics of generic au- 
tonomous quantum many-body systems on the energy 
shell in a similar way as 1-dim quantum maps describe 
(quantum) Poincare sections of 2-dim quantum (chaotic) 
systems (see e.g. B). 
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The total number of particles N = J^j n j 1S conserved, 
so the map U acts over Hilbert (Fock) space Ti. of di- 
mension TV — („) . The dynamics of a given initial 
many body state IV'(O)), which is an iteration of the map 
\ip(m)) = U\ip(m — 1)) = U m \ip(0)), can be performed 
most efficiently by observing that the kinetic part Ut is 
diagonal in the momentum basis \k) = c k ■ ■ ■ c£ |0), k\ < 
... < kpi while the potential part Uy is diagonal in 
the position basis \j) = c' j± ■ ■■ c] N \0},ji < ... < j N . 
The transformation between the two, Fj^ = (j\k), is 
an antisymmetrized N-dim discrete Fourier transforma- 
tion (DFT) on L-sites which has been efficiently coded 
in ~ TVlog 2 TV floating point operations (FPO) by fac- 
torizing L— site DFT to the product of O (L log 2 L) 2- 
site transformations parametrized with 2x2 sub-matrices 
(a,/3;j,S)jf, which are successively applied to creation 
operators, (ct , cl ) <— (act + /Jet, , jcj + 6c'-, ) , in all slater 

determinants n„cj n |0) which contain a particle at sites 
j or j'. Our algorithm (fermionic FFT) requires almost 
no extra storage apart from a vector of TV c-numbers and 
works for lattices of sizes L = 2P, 10, 12, 15, 20, 24, 30, 40. 
Therefore, the map (|^) is iterated on a vector tp^(m) — 

(k\ip(m)), using the matrix composition U = F*UvFUt 
in roughly 27V log 2 TV FPO per time step which is by far 
superior to complete diagonalization techniques (O (A/" 3 ) 
FPO), even for long time scales m = O (A/") when quan- 
tum dynamics becomes quasiperiodic due to discreteness 
of the spectrum of U. 
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FIG. 1. Current autocorrelation function Cj(m) against 
discrete time m for quantum ergodic (t = V = 4, lower set 
of curves for various sizes L) and mixed regime (t = V = 1, 
upper set of curves) with density p = ~. Averaging over 
entire Fock space is performed, M' = N , for L < 20, whereas 
random samples of J\f' = 12000, and J\f' — 160 initial states 
have been used for L = 24, and L = 32, respectively. 

First we consider decay of current autocorrelation 
function Cj(m) = (l/L)(J(m) J(0)) where J(m) = 
U fm JU m , and (.) = (l/N)Tr(.) is a 'microcanonical av- 
erage'. Note that J is diagonal in the momentum basis 
J\k) = J%\k), and (J) = 0. So Cj(m) can be evaluated 



by means of time-evolution of momentum initial states 



Cj(m) = LTV 5- J *'5Z J kPkk'( m ) 



(3) 



where pfe(m) = |(fc|^(m))| 2 = \{k\U m \k')\ 2 . For large 
sizes L, a smaller but uniformly random sample of TV"' 
initial states \k'), 1 <C TV' -C TV, is used in order to save 
computer time. Direct computation of Cj(m) for m < M 
can be performed in ~ (2MTVTV'/L) log 2 TV FPO, since 
due to translational symmetry one can simultaneously 
simulate the dynamics of L different states with different 
values of the conserved total momentum K = ^ n k' n 
(mod L). Using the eigenphases rj n and eigenstates \n) 
of evolution operator U, U\n) — e~' inn \n),n = 1...TV, 
one can write dissipative dc conductivity of such a kicked 



system a := J2n=i( d 4>Vn) 2 
Note that d<pr] n = (n\J\n). 
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FIG. 2. Stiffness Dj vs. 1/L at constant density p = ~ and 
for different values of control parameters in ergodic, t = V — 4 
and t — V = 2, and mixed, t = 1, V = 2 and t = V = 1, 
regime. Other parameters are the same as in fig.l. 

In fig.l we present numerical computation of corre- 
lation function Cj(m) for parameters t = V = 1 and 
t = V = 4, for various sizes L, but at fixed density 
p = N/L = j. Quite generally, Cj(m) exhibits fast 
relaxation on a time scale M* which is typically small, 
M * ~ 10, and roughly independent of L, and afterwards 
it fluctuates around averaged limiting value, the stiffness 



Dj = lim 

M 



1 

X. M 



M 



(4) 



m — 1 



where the strength of fluctuations decreases with in- 
cresing size L. Note again that TL L — > oo should be 

taken prior to the time-limit, lirnM->oo(l/-&0 X)m=i(-); 
which is for systems of finite size L here and below es- 
timated numerically as (1/M') X^m=M'+i(0 w i tn suffi- 
ciently large but fixed averaging time M' > M*; we 
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take M' = 30. If the system is quantum ergodic (case 
t = V = 4 of fig . 1 ) , Dj goes to zero and a remains finite 
as L — > oo (TV — > oo) and p = N/L fixed, whereas in the 
other case (t = V = 1 of fig.l) Dj remains well above 
zero as we approach TL whereas conductivity a diverges 
0. In fig. 2 we have analyzed the scaling of Dj with 
1/L. For large values of parameters, say t = V = 4, Dj 
is practically zero already for L s=s 20, while for smaller 
(but not small) control parmeters (t, V), Dj w Dy+f3/L 
where Df > 0. In close-to-critical case t = V = 2, we 
find larger correlation time M* ~ 10 2 , and hence use 
longer averaging time M' = 200. In fig. 3 we illustrate 
an ideal transport for t ~ V~ ~ 1 by plotting a persistent 
current M = Umjw_ >00 (l/M) X)m=i(^'l^( m )l^') vs - tne 
initial current Jg, . The normal transport in the ergodic 
regime £ = V = 4 is characterized by J? =0, while for 
t ~ V ~ 1 we find ideal transport with the persistent cur- 
rent being proportional to the initial current, JS = a Jg,. 
Proportionality constant a can be computed from (jj) 
Dj = (l/L)(j£,Jf,) = (a/L){,P), so a = 2Dj/[p(l- p)}, 
where (J 2 ) is given below (|J). 



Again similar behavior is found numerically for R(m) as 
for Cj(m) (the two quantities can be computed simulta- 
neously at no extra cost), namely it typically saturates 
within the same (short) correlation time M* to a roughly 

constant value R = limjvf->oo(l/-&0 Ylm=i R( m ) 
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FIG. 4. Relative localization dimension in Fock space, 
R(m) for data of fig.l. 




FIG. 3. Persistent current JS against initial current Jg (av- 
eraged over bins of size A J = 0.005) in the ergodic t — V = 4, 
(nearly) ergodic t = V = 2, and mixed regime, t = V = 1 and 
t = 1, V = 2. In all cases, L — 24 and p = i. 

Due to translational symmetry the total momentum K = 
y\ fcrife (mod L) is the only conserved quantity (appart 
from AT and parity), so the evolution of initial momentum 
state \k') takes place in Af K » AA/L dim. subspace Wjf, 
spanned by with K = \k\ := '}2, n k n . Starting with 
a momentum state fc'), the number of 'excited' states 
|fc) after time m is characterized by information entropy 

H (see also f|) as exp(- Y,%P%%>{™) lo gPj?fe'( TO ))- Aver- 
aging the entropy over uniformly random sample of A/*' 

initial states \k') we define relative localization dimension 
in Fock space as a measure of quantum ergodicity. 



R{m) = jj. exp 



I ^2 ^2 p kk> ( 

\ k> k 




m) log p£p(m) 
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FIG. 5. Limiting rel.loc.dim. R vs. 1/L for data of fig. 2. 

If there are no additional, even approximate conser- 
vation laws the blocks U m \u K may be modeled by cir- 
cular orthogonal ensemble (COE) of random matrices 
for sufficiently large m giving the maximal asymptotic 
(as Af — > oo) value of relative localization dimension, 
Rcoe ~ 0.655. This case corresponds to quantum er- 
godicity since p^, (m) , for m > M* , become pseudo ran- 
dom and independent of k and k' , hence the correlation 
function @ factorizes and yields C(m) = (J) 2 = 0. In- 
deed, as we show in fig. 4, such behavior is obtained only 
for sufficiently large parameters, say t = V = 4, while 
for smaller values of parameters t, V, R{m) saturates to 
a smaller value indicating that there may exist approx- 
imate conservation laws causing nontrivial localization 
inside the Fock space. Scaling with 1/L suggests that 
even TL of R is smaller than Rcoe for t ~ V ~ 1 (fig. 5). 
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FIG. 6. Steady-state current distribution divided by a 
gaussian P{I, Io)/P g auss{I) averaged over 279 initial orbits 
with |Jo| < 0.08 in the ergodic, t = V = 4, and mixed, 
t = V = 1, regime, and the finite-size microcanonical cur- 
rent distribution P mc (I). L = 24, p = |. 

Finally, we discuss current fluctuations, or more gener- 
aly, current distribution P^,(I) = (ip\S(I — J)\ip) giving a 
probability density of having a current / in a state \tp). 
We let the state tp with a 'good' known initial current 
lo to evolve for a long time from which we compute a 
stady-state current distribution (SSCD) 

M 

P(I; I ) = Jim - <*W> J(W(I Am))). 

Of course, delta functions should have a finite small 
width providing averaging over several states with 
J^k, I . In the quantum ergodic regime all states even- 
tually become populated, so SSCD P(I;Io) should be 
independent of initial current Iq and equal to the mi- 
crocanonical current distribution P mc {T) = (5(1 — J)). 
It has been shown by elementary calculation that in TL 
the latter becomes a gaussian, P mc (I) — > P ga uss(I) = 
{1/ ^2tt(J 2 )) cxv(-\l 2 /(J 2 )), while at any finite size L 
the first few moments are: 



(J 2 ) 



N(L - N) 



y(i- P )L, 



2(L-1) 

(J 4 ) _ 3(L- l)(2iV(L- N) -L) 
(J2)2 - 2N(L-2)(L-N) 
3(2p(l - p) - 1) 1 
2p(l -p) L 



(5) 



3 + 



+ 0(L- 2 ). 



Numerical results for L = 24 (see fig. 6) indicate that 
in the ergodic regime, t = V = 4, SSCD is al- 
ready in good agreement with microcanonical distribu- 
tion P mc (I), while in non-ergodic (mixed) regime, t = 
V = 1, SSCD is localized on a smaller range indicating 
that the current fluctuation is smaller than (J 2 ). Note 
that the mean I = J dIIP(I; Iq) is just a persistent cur- 
rent, so I = alo (see fig. 3). 



In this Letter we have presented numerical evidence, 
based on efficiently coded time evolution of a kicked 
fermionic system, in support of hypothesis that mixed 
(neither integrable nor ergodic) behavior of a quantum 
many-body system may survive TL provided that control 
parameters are not too far away from integrable points. 
It has been shown that in this regime ideal transport is 
possible. However, if the control parameters are suffi- 
ciently large we recover quantum ergodycity compatible 
with random matrix theory and normal transport prop- 
erties. It is interesting to note that at the transition 
point between the two regimes, where order parameter 
- stiffness Dj\l =oc (inferred from 1/L scaling) touches 
zero, the correlation time scale M* drastically increases 
what is reminiscent of a kind of dynamical phase transi- 
tion. Although only data for quarter-filled lattice (p = j) 
are presented here, we should stress that essentially the 
same conclusions follow from our data for other densities, 
p = !)§)§)§) with a general rule, that the border of 
quantum ergodic regime moves to slightly smaller values 
of control parameters t, V as the density p approaches 
|. It should be noted that statistics of eigenphases of 
evolution operator U has been computed as well and it 
has been found that in the ergodic regime level statis- 
tics is indeed that of COE while in the mixed regime it 
smoothly interpolates between Poisson and COE. 

Discussions with Prof. P. Prelovsck, and the financial 
support by Ministry of Science and Technology of the 
Republic of Slovenia are gratefully acknowledged. 
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